This R code includes the analyses and figures of the SMART manuscript based on the validation cohort of 95 AML patient samples.
It includes analyses of Junyan Lu, with extensions and modifications by Peter-Martin Bruch.
Load libraries and data
Libraries
library(table1)
##
## Attache Paket: 'table1'
## Die folgenden Objekte sind maskiert von 'package:base':
##
## units, units<-
library(pheatmap)
library(RColorBrewer)
library(ggbeeswarm)
## Lade nötiges Paket: ggplot2
## Warning: Paket 'ggplot2' wurde unter R Version 4.2.3 erstellt
library(grImport)
## Warning: Paket 'grImport' wurde unter R Version 4.2.3 erstellt
## Lade nötiges Paket: grid
## Lade nötiges Paket: XML
## Warning: Paket 'XML' wurde unter R Version 4.2.3 erstellt
library(cowplot)
library(geometry)
## Warning: Paket 'geometry' wurde unter R Version 4.2.3 erstellt
library(tidyverse)
## Warning: Paket 'tidyverse' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'tibble' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'readr' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'dplyr' wurde unter R Version 4.2.3 erstellt
## Warning: Paket 'lubridate' wurde unter R Version 4.2.3 erstellt
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ purrr 1.0.1 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ lubridate::stamp() masks cowplot::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Load processed data
load("../../data/Validation_data.Rdata")
filepath<-"./00_raw_data/"
source("./functions.R")
## Lade nötiges Paket: ggpubr
## Warning: Paket 'ggpubr' wurde unter R Version 4.2.3 erstellt
##
## Attache Paket: 'ggpubr'
## Das folgende Objekt ist maskiert 'package:cowplot':
##
## get_legend
##
## Attache Paket: 'survminer'
## Das folgende Objekt ist maskiert 'package:survival':
##
## myeloma
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
fontsize=18
Functions for patient characteristics table
my.render.cont <- function(x) {
with(stats.apply.rounding(stats.default(x), digits=3), c("",
"Median [Min, Max]"=sprintf("%s [%s, %s]", MEDIAN, MIN, MAX)))
}
my.render.cont_IQR <- function(x) {
with(stats.apply.rounding(stats.default(x), digits=3), c("",
"Median [Min, Max]"=sprintf("%s [%s, %s]", MEDIAN, IQR)))
}
my.render.cat <- function(x) {
c("", sapply(stats.default(x), function(y) with(
y,sprintf("%d (%0.0f %%)", FREQ, PCT))))
}
Preprocessing
Cytarabin screened twice with the concentration of 20000, average the values
avgTab <- group_by(screenData, name, Concentration1, Concentration2, patID) %>%
summarise(normVal = mean(normVal))
## `summarise()` has grouped output by 'name', 'Concentration1', 'Concentration2'.
## You can override using the `.groups` argument.
screenData <- distinct(screenData, Drug1, Concentration1, Drug2, Concentration2, name, wellType, patID) %>%
left_join(avgTab, by = c("patID","name","Concentration1","Concentration2"))
Calculate AUC for single drugs
Function for calculating AUC using trapezoidal rule
calcAUC <- function(value, conc) {
#value is a vector of normalized viability values
#conc is a vector of concentrations.
valueConc <- data.frame(viab=value, conc=conc)
valueConc <- valueConc[order(valueConc$conc),]
areaTotal <- 0
for (i in seq(length(value)-1)) {
areaEach <- (valueConc$viab[i] + valueConc$viab[i+1])*log(valueConc$conc[i+1]/valueConc$conc[i])*0.5
areaTotal <- areaTotal + areaEach
}
areaTotal/log(valueConc[nrow(valueConc),]$conc/valueConc[1,]$conc)
}
Calculate AUC
viabTab <- screenData %>%
filter(wellType == "single") %>%
mutate(Drug = Drug1,
conc = Concentration1)
aucTab.single <- group_by(viabTab, patID ,Drug) %>%
summarise(auc = calcAUC(normVal, conc))
## `summarise()` has grouped output by 'patID'. You can override using the
## `.groups` argument.
Calculate volumn under surface for combinations
getVolume=function(df) {
#find triangular tesselation of (x,y) grid
res=delaunayn(as.matrix(df[,-3]), output.options = TRUE)
#calulates sum of truncated prism volumes
sum(mapply(function(triPoints,A) A/3*sum(df[triPoints,3]),
split.data.frame(res$tri,seq_along(res$areas)),
res$areas))
}
calcVUS <- function(conc1, conc2, normVal) {
eachTab <- tibble(conc1 = conc1, conc2=conc2, normVal=normVal)
df <- eachTab %>%
mutate(conc1=log10(conc1), conc2 = log10(conc2)) %>%
mutate(conc1 = ifelse(is.infinite(conc1),0,conc1),
conc2 = ifelse(is.infinite(conc2), 0, conc2)) %>%
data.frame()
dfMax <- df
dfMax$normVal <- 1
v1 <- getVolume(df)
v2 <- getVolume(dfMax)
return(v1/v2)
}
#get combination
viabTab.combi <- screenData %>%
filter(wellType == "combi")
#calculate volunm under surface
aucTab.combi <- group_by(viabTab.combi, name, patID) %>%
summarise(auc = calcVUS(Concentration1, Concentration2, normVal)) %>%
dplyr::rename(Drug = name)
## `summarise()` has grouped output by 'name'. You can override using the
## `.groups` argument.
Combine auc of single and combi drugs to one table
aucTab <- bind_rows(aucTab.single, aucTab.combi)
Individual concentrations for visualization
concTab <- screenData %>%
filter(wellType %in% c("single","combi"))
Annotate in-vivo response and eln risk
aucTab <- left_join(aucTab, patMeta, by = "patID")
concTab <- left_join(concTab, patMeta, by = "patID")
Overview of the dataset
Pie-chart for ELN risk and in-vivo response
library(webr)
## Warning: Paket 'webr' wurde unter R Version 4.2.3 erstellt
plotTab <- aucTab %>% distinct(patID, ELN22, group) %>%
group_by(ELN22) %>%
mutate(n1 = length(patID)) %>%
mutate(ELNrisk = sprintf("%s\n(n=%s)",ELN22,n1)) %>%
group_by(group, ELNrisk) %>%
summarise(n = length(patID)) %>%
ungroup() %>%
dplyr::rename(invivo = group)
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
PieDonutCustom(plotTab, aes(ELNrisk, invivo, count=n, pieLabel = "ELN risk"),
showPieName = FALSE,
showRatioDonut = TRUE, showRatioPie = FALSE,
r0 = getOption("PieDonut.r0", 0.2),
r1 = getOption("PieDonut.r1", 0.7),
r2 = getOption("PieDonut.r2", 1), piLabel = "ELN risk", donutLabel = "in vivo response",
donutLabelSize = 4)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
aucTab %>% distinct(patID, ELN22, group) %>%
group_by(ELN22, group) %>%
count()
## # A tibble: 6 × 3
## # Groups: ELN22, group [6]
## ELN22 group n
## <chr> <fct> <int>
## 1 adverse Non-Responder 28
## 2 adverse Responder 8
## 3 favorable Non-Responder 6
## 4 favorable Responder 10
## 5 intermediate Non-Responder 14
## 6 intermediate Responder 29
writexl::write_xlsx(rownames_to_column(as.data.frame(plotTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4A",".xlsx"))
Patient characteristics table
patMeta_Validation<-left_join(patMeta_Validation, patMeta, by="patID")
label(patMeta_Validation$AGE) <-"Age (years)"
label(patMeta_Validation$diagnosis) <-"Diagnosis"
label(patMeta_Validation$trtstr) <-"Prescribed treatment"
label(patMeta_Validation$SEX) <-"Sex"
label(patMeta_Validation$BMB) <-"Tumor cell infiltration"
label(patMeta_Validation$group) <-"In vivo response group"
label(patMeta_Validation$ELN22) <-"ELN-22 risk group"
# table 1 printing
table1(~ SEX + AGE + diagnosis + BMB + trtstr + group + ELN22, data=patMeta_Validation, overall="Total", render.continuous=my.render.cont, render.categorical=my.render.cat)
Supp_Table_Val_Cohort<-patMeta_Validation %>%
rename(PatientID=patID, Age=AGE, Sex=SEX, Treatment=trtstr, `Tumor cell infiltration`=BMB, Diagnosis=diagnosis, `In vivo response group`=group, `ELN-22 risk group`=ELN22 ) %>%
select(-EFSSTAT, -EFSTM, -EFSTMU) %>%
select(PatientID, Age, Sex, Treatment, `Tumor cell infiltration`, Diagnosis, `In vivo response group`, `ELN-22 risk group`)
writexl::write_xlsx(Supp_Table_Val_Cohort,path="../../docs/table_S3_Validation_patients.xlsx")
Drug response heatmap
colAnno <- aucTab %>% distinct(patID, ELN22, group) %>%
mutate(ELN22 = factor(ELN22, levels = c("adverse","intermediate","favorable")),
group = factor(group, levels = c("Non-Responder","Responder"))) %>%
arrange(ELN22, group) %>%
dplyr::rename(`ELN risk` = ELN22, `in vivo response`= group) %>%
column_to_rownames("patID")
annoColor <- list(`ELN risk` = c(adverse = "#E18727FF", intermediate = "#6F99ADFF", favorable = "#20854EFF"),
`in vivo response` = c(`Non-Responder` = "orangered3", Responder = "deepskyblue4"))
viabMat <- viabMat[,rownames(colAnno)]
pheatmap(viabMat, annotation_col = colAnno,
cluster_cols = FALSE, annotation_colors = annoColor,
# filename=paste0("../../docs/validation_drugHeatmap.pdf"),
show_colnames = FALSE,
treeheight_row = 20)
With hierarchical clustering
viabMat <- viabMat[,rownames(colAnno)]
viabMat.scaled <- mscale(viabMat, censor = 4)
p <- pheatmap(viabMat.scaled, annotation_col = colAnno,
cluster_cols = TRUE, annotation_colors = annoColor,
color=colorRampPalette(c("deepskyblue3", "white", "orangered4"))(100), # colorRampPalette(c(blue2, "white", red2))(100)
clustering_method = "ward.D2",
main = "Hierarchical clustering of ex vivo responses in AML validation cohort",
show_colnames = FALSE,
treeheight_row = 20, silent = TRUE)
cowplot::plot_grid(p$gtable)
# ggsave("../../docs/HC_heatmap_Validation.png", plot = p$gtable, height = 8, width =15)
ggsave("../../docs/Figure_ED4.pdf", plot = p$gtable, height = 8, width =15)
writexl::write_xlsx(rownames_to_column(as.data.frame(viabMat.scaled),var="rowname"), path=paste0("../../","docs/Source data/","ED4",".xlsx"))
Test correlations between ex-vivo response and in-vivo responses
Univariate test (without blocking for ELN risk)
testTab <- aucTab
resTab.uni <- group_by(testTab, Drug) %>% nest() %>%
mutate(m = map(data, ~t.test(auc~group,., var.equal=TRUE))) %>%
mutate(res = map(m, broom::tidy)) %>%
unnest(res) %>%
select(Drug, estimate, p.value) %>%
arrange(p.value) %>% ungroup() %>%
mutate(p.adj = p.adjust(p.value, method ="BH"))
sigDrugs <- filter(resTab.uni, p.adj <= 0.1)$Drug
Volcano plot
fdrCut=0.1
plotTab <- resTab.uni %>%
mutate(dir = ifelse(estimate > 0, "more resistant", "more sensitive")) %>%
mutate(dir = ifelse(p.adj > fdrCut, "n.s.", dir))
pCut <- -log10(min(filter(resTab.uni, p.adj > fdrCut)$p.value))
pVolcano <- ggplot(plotTab, aes(x=estimate, y=-log10(p.value))) +
geom_point(aes(fill = dir), shape = 21) +
scale_fill_manual(values = c(`more resistant` = "orangered3", `more sensitive` = "deepskyblue3", n.s. = "grey80"), name = "") +
geom_hline(yintercept = pCut, color = "grey50", linetype = "dashed") +
ggrepel::geom_text_repel(data = filter(plotTab, p.adj <= fdrCut),aes(label = gsub("_"," + ",Drug)), size=5) +
ylab(expression(-log[10]*'('*italic(P)~value*')')) + xlab("Mean viability difference") +
ggtitle(bquote(italic("In vivo")~"response"~"~"~italic("ex vivo")~"response")) +
theme_classic()+
theme(plot.title = element_text(hjust=0.5, face = "bold", size=18),
axis.title = element_text(size=16),
axis.text = element_text(size=16),
legend.text = element_text(size=14),
legend.background = element_rect(fill = NA),
axis.line = element_blank(),
legend.position = c(0.75,0.2),
plot.margin = margin(0,20,0,20),
panel.border = element_rect(size=1.5, fill = NA)
)
pVolcano
Boxplot
resTab.sig <- filter(resTab.uni, p.adj <= 0.1)
plotList_box <- lapply(seq(nrow(resTab.sig)) ,function(i){
rec <- resTab.sig[i,]
#get the df
df_plot <- filter(aucTab, Drug == rec$Drug)
pVal <- formatC(rec$p.value, digits = 2)
drugName <- rec$Drug
#make plot
p <- ggplot(df_plot, aes(x=group, y=auc)) +
geom_boxplot(aes(fill = group), width = 0.8, size=0.8, alpha = 0.85,
color = "black", outlier.shape = NA) +
scale_fill_manual(values = c("Non-Responder"="orangered3","Responder"="deepskyblue4")) +
geom_beeswarm(cex=3, dodge.width=0.8, shape=21,
fill="slategray4", size=4, alpha= 0.7) +
geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
coord_cartesian(ylim = c(0,1.5)) +
#ggtitle(sprintf("%s (P = %s)",rec$Drug, ) +
ggtitle(bquote(.(drugName)~"("*italic("P")~"="~.(pVal)~")")) +
#xlab(expression(paste(italic("in vivo")," response"))) +
xlab("") +
ylab("Viability") +
guides(color="none", fill="none") +
theme_classic()+
theme(plot.title = element_text(hjust=0, face = "bold", size=18),
axis.title = element_text(size=16),
axis.text = element_text(size=16),
legend.text = element_text(size=14),
legend.background = element_rect(fill = NA),
axis.line = element_blank(),
legend.position = c(0.75,0.2),
plot.margin = margin(0,20,0,20),
panel.border = element_rect(size=1.5, fill = NA)
)
return( p )
})
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
names(plotList_box) <- resTab.sig$Drug
pBox <- cowplot::plot_grid(plotlist = plotList_box, ncol=2)
pBox
writexl::write_xlsx(rownames_to_column(as.data.frame(aucTab),var="rowname"), path=paste0("../../","docs/Source data/","ED8A",".xlsx"))
Dose-response curves for selected drugs
example_drugs <- sigDrugs[sigDrugs!= "Daunorubicin_Cytarabin"]
plotList_DResp <- lapply(example_drugs ,function(drugname){
plotTab <- filter(concTab, Drug1 == drugname, Drug2 == "DMSO")
ciTab <- group_by(plotTab, group, Concentration1) %>%
summarise(sdVal = sd(normVal), n=length(patID), normVal = mean(normVal)) %>%
mutate(se = sdVal/sqrt(n), ci = 1.96*se)
#plot values: real with dots, predicted with line
ggplot(plotTab,aes(x= Concentration1, y=normVal, color=group, fill=group)) +
geom_point(shape = 21, size=3.5, alpha = 0.2, color="black") +
geom_errorbar(data = ciTab, aes(x=Concentration1, y=normVal, ymax = normVal + ci, ymin=normVal-ci), width=0.1, size=1) +
scale_color_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered3")) +
scale_fill_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered3")) +
coord_cartesian(ylim = c(0,1.5)) +
geom_smooth(data = ciTab, method = "fitIC50", se = FALSE) +
scale_x_continuous(trans='log10') +
labs(title = paste0("",drugname) )+
xlab(expression("Concentration"~(mu*M))) + ylab("Viability") +
guides(fill="none", color="none") +
theme_classic()+
theme(plot.title = element_text(hjust=0, face = "bold", size=18),
axis.title = element_text(size=16),
axis.text = element_text(size=16),
legend.text = element_text(size=14),
legend.background = element_rect(fill = NA),
axis.line = element_blank(),
legend.position = c(0.75,0.2),
plot.margin = margin(0,20,0,20),
panel.border = element_rect(size=1.5, fill = NA)
)
})
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'group'. You can override using the
## `.groups` argument.
pDose <- cowplot::plot_grid(plotlist = plotList_DResp, ncol=2)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
pDose
writexl::write_xlsx(rownames_to_column(as.data.frame(concTab),var="rowname"), path=paste0("../../","docs/Source data/","ED8B",".xlsx"))
Multi-variate test (with blocking for ELN risk)
testTab <- filter(aucTab, Drug %in% sigDrugs)
resTab.multi <- group_by(testTab, Drug) %>% nest() %>%
mutate(m = map(data, ~car::Anova(lm(auc~group + ELN22,.)))) %>%
mutate(res = map(m, broom::tidy)) %>%
unnest(res) %>% filter(term == "group") %>%
select(Drug, p.value) %>%
arrange(p.value) %>% ungroup() %>%
mutate(p.adj = p.adjust(p.value, method ="BH"))
head(resTab.multi)
## # A tibble: 6 × 3
## Drug p.value p.adj
## <chr> <dbl> <dbl>
## 1 Vindesine 0.000945 0.00330
## 2 Fludarabine 0.00152 0.00330
## 3 Vincristine 0.00163 0.00330
## 4 Cladribine 0.00189 0.00330
## 5 NSC 652287 0.0133 0.0186
## 6 Daunorubicin_Cytarabin 0.0465 0.0542
Test associations between ex-vivo and in-vivo responses, within each ELN group
testTab <- filter(aucTab, Drug %in% resTab.sig$Drug)
resTab <- testTab %>%
group_by(Drug, ELN22) %>% nest() %>%
mutate(m = map(data, ~t.test(auc~group,.,var.equal=TRUE))) %>%
mutate(res = map(m, broom::tidy)) %>%
unnest(res) %>%
select(Drug, p.value, estimate) %>%
arrange(p.value) %>% ungroup() %>%
mutate(p.adj = p.adjust(p.value, method ="BH"))
## Adding missing grouping variables: `ELN22`
head(resTab)
## # A tibble: 6 × 5
## ELN22 Drug p.value estimate p.adj
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 adverse Vindesine 0.00331 0.249 0.0331
## 2 adverse Fludarabine 0.00417 0.136 0.0331
## 3 adverse Vincristine 0.00473 0.208 0.0331
## 4 adverse Cladribine 0.00663 0.233 0.0348
## 5 adverse Daunorubicin_Cytarabin 0.0138 0.190 0.0579
## 6 adverse NSC 652287 0.0443 0.183 0.155
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","ED9B",".xlsx"))
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4C",".xlsx"))
There are significant p-values, but p.adj is borderline significant
Plot p-value heatmap
plotTab <- resTab %>% mutate(pSign = -log10(p.value) * sign(estimate)) %>%
mutate(ifSig = case_when(p.adj <=0.1 ~ "*",
TRUE ~ ""))
pPheatmap <- ggplot(plotTab, aes(x=factor(Drug), y=factor(ELN22, levels = c("adverse", "intermediate", "favorable")), fill = pSign)) +
geom_tile() +
geom_text(aes(label = ifSig), col = "black", size=8, face="bold") +
scale_fill_gradient2(low = "blue", mid = "white",high = "red", name = bquote("-log10("*italic("P")*")")) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
xlab("") + ylab("ELN risk group") +
theme(axis.title = element_text(size=18),
axis.text = element_text(size=16),
plot.margin = margin(5,20,5,20))
## Warning in geom_text(aes(label = ifSig), col = "black", size = 8, face =
## "bold"): Ignoring unknown parameters: `face`
pPheatmap
Boxplot of significant associations
resTab.sig <- filter(resTab, p.adj <= 0.1)
pList <- lapply(unique(resTab.sig$Drug), function(dd) {
plotTab <- filter(testTab, Drug == dd) %>%
mutate(ELN22 = factor(ELN22, levels = rev(c("favorable","intermediate","adverse"))))
pTab <- filter(resTab, Drug == dd) %>%
mutate(pText = sprintf("italic(P)~'='~%s",formatC(p.value, digits=2)))
ggplot(plotTab, aes(x=group, y=auc)) +
geom_boxplot(aes(fill = group), width = 0.8, size=0.8, alpha = 0.85, color = "black", outlier.shape = NA) +
scale_fill_manual(values = c("Non-Responder"="orangered3","Responder"="deepskyblue4")) +
geom_beeswarm(cex=3, dodge.width=0.8, shape=21,
fill="slategray4", size=4, alpha= 0.7) +
geom_text(data = pTab, aes(label = pText), x=2.1,y=1.3, vjust=1, parse= TRUE,size=5)+
geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
coord_cartesian(ylim = c(0,1.5)) +
ggtitle(gsub("_"," + ",dd)) +
xlab(expression(paste(italic("in vivo")," response"))) +
ylab("Viability (AUC)") +
guides(color="none", fill="none") +
theme_classic()+
facet_wrap(~factor(ELN22, levels = rev(c("favorable","intermediate","adverse")))) +
theme(plot.title = element_text(hjust=0, face = "bold", size=18),
axis.title = element_text(size=16),
axis.text = element_text(size=16),
axis.text.x =element_text(angle = 45, hjust = 1),
legend.text = element_text(size=14),
legend.background = element_rect(fill = NA),
axis.line = element_blank(),
legend.position = c(0.75,0.2),
plot.margin = margin(0,20,0,20),
panel.border = element_rect(size=1.5, fill = NA),
strip.text.x = element_text(size = 16)
)
})
names(pList) <- unique(resTab.sig$Drug)
#cowplot::plot_grid(plotlist = pList, ncol=1)
mainDrug <- c("Vincristine","Daunorubicin + Cytarabin", "Daunorubicin_Cytarabin")
pList$Daunorubicin_Cytarabin$labels$y<-"Viability (VUC)"
pList$Daunorubicin_Cytarabin$labels$title<-"Daunorubicin + Cytarabine"
pBoxELN_main <- plot_grid(plotlist = pList[names(pList) %in% mainDrug],ncol=1)
pBoxELN_supp <- plot_grid(plotlist= pList[!names(pList) %in% mainDrug],ncol=1)
pBoxELN_main
Dose-response curves for selected drugs per elnGroup
example_drugs <- sigDrugs[sigDrugs!= "Daunorubicin_Cytarabin"]
plotList_DResp <- lapply(example_drugs ,function(drugname){
plotTab <- filter(concTab, Drug1 == drugname, Drug2 == "DMSO") %>%
mutate(ELN22 = factor(ELN22, levels=c("adverse","intermediate","favorable")))
meanTab <- group_by(plotTab, Drug1, ELN22, group, Concentration1) %>%
summarise(normVal = mean(normVal, na.rm=TRUE))
#plot values: real with dots, predicted with line
ggplot(plotTab,aes(x= Concentration1, y=normVal, color=group, fill=group)) +
geom_point(shape = 21, size=3.5, alpha = 0.5, color="black") +
scale_color_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered4")) +
scale_fill_manual(values = c("Responder" = "deepskyblue4", "Non-Responder"="orangered4")) +
coord_cartesian(ylim = c(0,1.5)) +
geom_smooth(data=meanTab, method = "fitIC50", se = FALSE) +
scale_x_continuous(trans='log10') +
labs(title = paste0("",drugname) )+
xlab(expression("Concentration"~(mu*M))) + ylab("Viability") +
facet_wrap(~ELN22) +
guides(fill="none", color="none") +
theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_blank(),
panel.border = element_rect(size=1.5, fill = NA),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2))
})
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Drug1', 'ELN22', 'group'. You can override
## using the `.groups` argument.
cowplot::plot_grid(plotlist = plotList_DResp, ncol=1)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Computation failed in `stat_smooth()`
## Caused by error in `loadNamespace()`:
## ! es gibt kein Paket namens 'dr4pl'
Associate drug responses with clinical outcomes
Prepare dataset
Outcome data
# clinicTab <- readxl::read_xlsx(paste0("../../data/Validation_cohort_survival.xlsx")) %>%
clinicTab <- patMeta_Validation %>%
select(patID, EFSSTAT, EFSTM) %>%
mutate(across(-patID, as.numeric)) %>%
dplyr::rename(patID = patID) %>%
mutate(patID = as.character(patID)) %>%
distinct(patID, .keep_all = TRUE)
#censor survival data at five years
clinicTab <- mutate(clinicTab, EFSSTAT = ifelse(EFSTM>=60,0,EFSSTAT), EFSTM = ifelse(EFSTM > 60, 60, EFSTM)) %>%
mutate(EFSTM = EFSTM/12)
ELN data
elnTab <- patMeta %>%
distinct(patID, ELN22) %>%
mutate(ELN22 = factor(ELN22, levels = c("favorable","intermediate","adverse")))
survTab <- left_join(elnTab, clinicTab, by = "patID")
Drug response data: use AUC for single drugs and individual concentrations for combinatorial drugs
drugTab <- aucTab %>%
filter(Drug %in% sigDrugs)
Test whether individual drug response associated with survival using univariate model
testTab <- survTab %>% pivot_longer(-c(patID,ELN22)) %>%
mutate(stat = ifelse(str_detect(name, "STAT"),"endpoint","time"),
outcome = str_remove(name, "STAT|TM")) %>%
select(-name) %>% filter(!is.na(value)) %>%
pivot_wider(names_from = stat, values_from = value) %>%
left_join(select(drugTab, patID, Drug, auc))
## Joining with `by = join_by(patID)`
resTab <- group_by(testTab, outcome, Drug) %>% nest() %>%
mutate(m = map(data, ~com(.$auc, .$time,.$endpoint))) %>%
unnest(m) %>%
arrange(p) %>%
select(outcome, Drug, p, HR, lower, higher) %>%
ungroup() %>%
mutate(p.adj = p.adjust(p, method ="BH"))
resTab.sig <- filter(resTab, p<0.05)
resTab.sig
## # A tibble: 5 × 7
## outcome Drug p HR lower higher p.adj
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 EFS Fludarabine 0.00728 13.7 2.03 93.1 0.0257
## 2 EFS Cladribine 0.0104 4.71 1.44 15.4 0.0257
## 3 EFS Vincristine 0.0167 6.75 1.41 32.3 0.0257
## 4 EFS Vindesine 0.0168 5.50 1.36 22.2 0.0257
## 5 EFS Daunorubicin_Cytarabin 0.0183 5.28 1.32 21.1 0.0257
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4D",".xlsx"))
pList <- lapply(seq(nrow(resTab.sig)), function(i) {
rec <- resTab.sig[i,]
plotTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug)
km(plotTab$auc, plotTab$time, plotTab$endpoint,
gsub("_"," + ",rec$Drug), pval = rec$p, stat = "maxstat", showTable = FALSE)
})
names(pList) <- unique(resTab.sig$Drug)
pKM_main <- pList[names(pList) %in% mainDrug]
pKM_supp <- pList[!names(pList) %in% mainDrug]
pList_withtable <- lapply(seq(nrow(resTab.sig)), function(i) {
rec <- resTab.sig[i,]
plotTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug)
km(plotTab$auc, plotTab$time, plotTab$endpoint,
rec$Drug, pval = rec$p, stat = "maxstat", showTable = TRUE)
})
names(pList_withtable) <- unique(resTab.sig$Drug)
pList_withtable[names(pList_withtable) %in% mainDrug]
## $Vincristine
##
## $Daunorubicin_Cytarabin
pList_withtable[!names(pList_withtable) %in% mainDrug]
## $Fludarabine
##
## $Cladribine
##
## $Vindesine
Supplementary Figure
cowplot::plot_grid(plotlist =pKM_supp, ncol=2)
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
Main figure panel
pLegend <- get_legend(pKM_main$Vincristine)
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
pKM_main$Daunorubicin_Cytarabin$labels$title<-"Daunorubicine + Cytarabine"
pKM_main <- plot_grid(plot_grid(pKM_main$Vincristine + theme(legend.position = "none"),
pKM_main$Daunorubicin_Cytarabin + theme(legend.position = "none"), ncol = 2),
pLegend, ncol=1 ,rel_heights = c(1,0.1))
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
Test whether individual drug response associated with survival using univariate model per ELN risk group
resTab <- group_by(testTab, outcome, Drug, ELN22) %>% nest() %>%
mutate(m = map(data, ~com(.$auc, .$time,.$endpoint))) %>%
unnest(m) %>%
arrange(p) %>%
select(outcome, Drug, p, HR, lower, higher) %>%
ungroup() %>%
mutate(p.adj = p.adjust(p, method ="BH"))
## Adding missing grouping variables: `ELN22`
resTab.sig <- filter(resTab, p<0.05)
resTab.sig
## # A tibble: 4 × 8
## ELN22 outcome Drug p HR lower higher p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 adverse EFS Cladribine 0.00355 17.1 2.54 115. 0.0541
## 2 adverse EFS Fludarabine 0.00515 86.9 3.81 1984. 0.0541
## 3 adverse EFS Vincristine 0.0229 14.4 1.45 143. 0.160
## 4 adverse EFS Daunorubicin_Cytarabin 0.0326 7.32 1.18 45.5 0.171
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","ED10",".xlsx"))
pList <- lapply(seq(nrow(resTab.sig)), function(i) {
rec <- resTab.sig[i,]
leg <- NULL
pELN <- lapply(c("adverse","intermediate","favorable"), function(eln) {
eachRec <- filter(resTab, ELN22 == eln, outcome == rec$outcome, Drug == rec$Drug)
plotTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug, ELN22== eln)
kmPlot <- km(plotTab$auc, plotTab$time, plotTab$endpoint,
paste0(eln," group"), pval = eachRec$p, stat = "maxstat", showTable = TRUE)
leg <<- cowplot::get_legend(kmPlot)
kmPlot <- kmPlot + theme(legend.position = "none")
kmPlot
})
gTitle <- ggplot() +
labs(title = rec$Drug) +
theme_void() + theme(plot.title = element_text(hjust = 0.5, face="bold", size= 15))
cowplot::plot_grid(gTitle,cowplot::plot_grid(plotlist = pELN, nrow=1), leg, ncol=1, rel_heights = c(0.1,1,0.2))
})
cowplot::plot_grid(plotlist = pList, ncol=1)
Test outcome prediction based on drug response plus ELN risk using multi-variate model
resTab <- group_by(testTab, outcome, Drug) %>% nest() %>%
mutate(m = map(data, ~coxph(
Surv(time, endpoint) ~
auc + ELN22,
data = .))) %>%
mutate(res=map(m, broom::tidy)) %>%
unnest(res) %>%
filter(term == "auc") %>%
arrange(p.value) %>%
select(outcome, Drug, p.value)
resTab.sig <- filter(resTab, p.value<0.05)
resTab.sig
## # A tibble: 5 × 3
## # Groups: outcome, Drug [5]
## outcome Drug p.value
## <chr> <chr> <dbl>
## 1 EFS Fludarabine 0.00537
## 2 EFS Cladribine 0.0185
## 3 EFS Vincristine 0.0211
## 4 EFS Daunorubicin_Cytarabin 0.0358
## 5 EFS Vindesine 0.0376
Plot significant associations with forest plots
plotHazard <- function(survRes, title = "") {
sumTab <- summary(survRes)$coefficients
confTab <- summary(survRes)$conf.int
#correct feature name
nameOri <- rownames(sumTab)
plotTab <- tibble(feature = c(rownames(sumTab), "reference (favorable)"),
HR = c(sumTab[,2],1),
p = c(sumTab[,5],1),
Upper = c(confTab[,4],1),
Lower = c(confTab[,3],1)) %>%
mutate(Upper = ifelse(Upper > 10, 10, Upper)) %>%
mutate(feature = case_when(feature == "auc" ~ "ex vivo resistance",
feature == "ELN22intermediate" ~ "intermediate",
feature == "ELN22adverse" ~ "adverse",
TRUE ~ feature)) %>%
mutate(pLabel = ifelse(!str_detect(feature, "reference"),
sprintf("italic(P)~'='~'%s'", formatNum(p, digits = 1)),
"")) %>%
mutate(feature = factor(feature, levels = c("reference (favorable)","intermediate","adverse","ex vivo resistance")))
ggplot(plotTab, aes(x=feature, y = HR)) +
geom_hline(yintercept = 1, linetype = "dotted", color = "grey50") +
geom_point(position = position_dodge(width=0.8), size=3, color = "black") +
geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.0, size=0.5,color = "grey20") +
geom_text(position = position_nudge(x = 0.3),
aes(y = HR, label = pLabel),
color = "black", size =5, parse = TRUE) +
expand_limits(y=c(-0.5,0))+
ggtitle(title) +
ylab("Hazard ratio") +
coord_flip() +
theme_full +
theme(legend.position = "none", axis.title.y = element_blank())
}
pList <- lapply(seq(nrow(resTab.sig)), function(i) {
rec <- resTab.sig[i,]
eachTab <- filter(testTab, outcome == rec$outcome, Drug == rec$Drug) %>%
select(time, endpoint, auc, ELN22) %>%
mutate(auc = (auc-mean(auc))/(2*sd(auc)))
survRes <- runCox(eachTab)
plotHazard(survRes, paste0(rec$outcome, " ~ ", rec$Drug))
})
names(pList) <- gsub("_"," + ", resTab.sig$Drug)
pList$`Daunorubicin + Cytarabin`$labels$title<-"EFS ~ Daunorubicin + Cytarabine"
writexl::write_xlsx(rownames_to_column(as.data.frame(testTab),var="rowname"), path=paste0("../../","docs/Source data/","Fig4E",".xlsx"))
pHR_main <- plot_grid(plotlist = pList[names(pList) %in% mainDrug], ncol=1)
pHR_supp <- pList[!names(pList) %in% mainDrug]
For Supplementary Figure
cowplot::plot_grid(plotlist= pHR_supp, ncol=2)
